#######################################
######## Joo et al. JCR Codes #########
#######################################


#############################
######### FIGURE 1 ##########
#############################
rm(list=ls())

library(foreign)
library(mvtnorm)
library(ggplot2)

wd <- setwd("Your_Directory/JCR_Jooetal_replication")


##### Figure 1a #####

dat<- read.dta("C:/Users/MinHyung_Joo/Dropbox/Research/Survey/Foreign Policy Project/Leader-level/JCR_replication_files/JCR_replication_March2024/JCR_Jooetal_replication/FP_data.dta")

rdat <- subset(dat, select=c('right2', 'year'))
rdat$count <- rep(1)
right_year <- aggregate(rdat, by=list(rdat$year), FUN=sum, rm.na=T)
right_year$right_pct <- as.vector(right_year[,2]/right_year[,4])*100

pdf(paste(wd, "/Figures/fig1a.pdf", sep=""), width=6,height=6, paper='special')

fig1a <- ggplot(right_year, aes(Group.1, right_pct)) + 
  geom_line(size=0.8) + labs(x ="Year", y = "Percentage of Right Wing Populist Leaders") +
  theme_bw()+
  geom_hline(yintercept=0, linetype="dashed") +
  theme(axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text = element_text(size = 12),
        panel.background = element_blank()) 
fig1a

dev.off()


##### Figure 1b #####

lcount <- read.dta("region_count.dta")

theTable <- as.data.frame(table(lcount$region8))

pdf(paste(wd, "/Figures/fig1b.pdf", sep=""), width=6,height=6, paper='special')

fig1b <- ggplot(theTable, aes(x=reorder(Var1, -Freq), Freq)) + 
  geom_bar(stat = "identity", width=0.5) + labs(x ="", y = "Number of Right Wing Populist Leaders") +
  theme_bw()+
  theme(axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text = element_text(size = 12),
        axis.text.x = element_text(angle = 45, hjust=1, face='bold'),
        panel.background = element_blank())

fig1b

dev.off()



######################
###### FIGURE 4 ######
######################
library(mvtnorm)
library(foreign)

set.seed(88110905)
coef <- c(1.698414,
          2.888438,
          -1.238306,
          -0.1168182,
          8.406813,
          1.27671,
          0.4112271,
          -0.2929019,
          0.1166458,
          0.0000485,
          0.184708,
          -0.0550502,
          0.3779452,
          0.0900058,
          -2.840332)

vcov <- as.matrix(read.csv("vcov_largeN.csv"))

dat <- read.dta("FP_data.dta")

lag <- 0
c <- median(dat$cinc, na.rm = T)
parti <- median(dat$partipdem, na.rm = T)
SD <- sd(dat$partipdem, na.rm=T)
r <- 0
maj <- 0
ally <- 1
nai <- median(dat$age, na.rm=T)
man <- 1
milcar <- 0
fight <- 0
milfight <- 0
edu <- 3
reb <- 0
irreg <- 0

B <- 1000

pardem <- seq(0,1, by =0.01) 

out1 <- matrix(NA, ncol=3, nrow=length(pardem))

for(i in 1:length(pardem)){ 
  pdem <- pardem[i] 	
  pr.bs <- matrix(NA,B,1)
  for(b in 1:B){ 
    newcoef <-  rmvnorm(1, mean=coef, sigma=vcov)
    bs.sum <- newcoef[1]*lag + newcoef[2]*0*pdem + newcoef[3]*pdem + newcoef[4]*0 + newcoef[5]*c +newcoef[6]*maj + newcoef[7]*ally
    + newcoef[8]*man + newcoef[9]*edu + newcoef[10]*milcar + newcoef[11]*fight + newcoef[12]*milfight + newcoef[13]*reb
    + newcoef[14]*irreg + newcoef[15]
    
    
    bs.sum2 <-  newcoef[1]*lag + newcoef[2]*1*pdem+ newcoef[3]*pdem + newcoef[4]*1 + newcoef[5]*c +newcoef[6]*maj + newcoef[7]*ally
    + newcoef[8]*man + newcoef[9]*edu + newcoef[10]*milcar + newcoef[11]*fight + newcoef[12]*milfight + newcoef[13]*reb
    + newcoef[14]*irreg + newcoef[15]
    pr.bs[b,] <- plogis(bs.sum2) - plogis(bs.sum) 
  }
  out1[i,] <- c(mean(pr.bs, na.rm=TRUE), quantile(pr.bs, c(0.025, 0.975), na.rm=TRUE)) 
  rm(bs.sum)
}

pdf(paste(wd, "/Figures/fig4.pdf", sep=""), width=6,height=6,paper='special')

par(mar=c(5,5,1.5,1.5))
plot(pardem, out1[,1], type="n", xlab="Participatory Democracy", ylab="Change in Pr(MIDS=1)",
     cex.lab=1.5, cex.axis=1.5, ylim = c(-0.3,0.8))
lines(lowess(pardem, out1[,1]), lty=1, lwd=2)
lines(lowess(pardem, out1[,2]), lty=2, lwd=2)
lines(lowess(pardem, out1[,3]), lty=2, lwd=2)
abline(h=0, lty=2, col='red', lwd=1)

dev.off()



# End of File
